Absorbing state phase transitions with quenched disorder 
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Quenched disorder - in the sense of the Harris criterion - is generally a relevant perturbation 
at an absorbing state phase transition point. Here using a strong disorder renormalization group 
framework and effective numerical methods we study the properties of random fixed points for sys- 
tems in the directed percolation universality class. For strong enough disorder the critical behavior 
is found to be controlled by a strong disorder fixed point, which is isomorph with the fixed point 
of random quantum Ising systems. In this fixed point dynamical correlations are logarithmically 
slow and the static critical exponents are conjecturedly exact for one-dimensional systems. The 
renormalization group scenario is confronted with numerical results on the random contact process 
in one and two dimensions and satisfactory agreement is found. For weaker disorder the numerical 
results indicate static critical exponents which vary with the strength of disorder, whereas the 
dynamical correlations are compatible with two possible scenarios. Either they follow a power-law 
decay with a varying dynamical exponent, like in random quantum systems, or the dynamical 
correlations are logarithmically slow even for weak disorder. For models in the parity conserving uni- 
versality class there is no strong disorder fixed point according to our renormalization group analysis. 



I. INTRODUCTION 

Nonequilibrium many-particle systems are often de- 
scribed by stochastic models in which some degrees of 
freedom with fast relaxation behavior are integrated out 
and are replaced by a noise term, which may represent 
thermal or quantum fluctuations, chaotic motion, etc. 
The classification of steady states of these nonequilib- 
rium models is an important task and we are currently 
witnessing considerable theoretical progress in this field. 
Of particular interest are systems that in the steady state 
show long-range spatial and temporal correlations, as is 
the case when they are in the vicinity of a nonequilibrium 
phase transition point. A special type of these are the 
transitions between an active phase and an inactive one 
where the particles are absorbed into a state without fluc- 
tuations. These absorbing state phase transitions play an 
important role in physics, chemistry and even biologji^. 

Classification of absorbing state phase transitions has 
shown that the universality class of directed percola- 
tion is particularly robust. It contains models with a 
scalar order parameter, absence of conservation laws, and 
short range interactions^^. Well known models with a 
phase transition in this universality class are the contact 
process^ and the Ziff-Gulari-Barshad model of catalytic 
reactions^. When there is a conservation law present, 
other universality classes can appear, the best known of 



which is the parity conserving class^ which includes the 
branching-annihilating random walk with even number 
of offspring^ and the generalized contact process with 
two different absorbing states^. Most recently reaction- 
diffusion models with multiple branching and annihila- 
tion processes were also intensively studiedSi. Sandpile 
models with conserved energy can also be related to ab- 
sorbing state phase transitions^''. 

Quenched, i.e. time-independent, disorder is an in- 
evitable feature of many real processes and could play 
an important role in stochastic particle systems too. As 
an example, it has been argued that due to the presence 
of some form of disorder the directed percolation univer- 
sality class has not yet been seen in real experimentsii, 
such as in catalytic reactions^, in dcpinning transitionsi^, 
and in the flow of granular matter^'^ (for a review, seeii). 
In stochastic particle systems, disorder is represented by 
position dependent reaction rates and its relevance can 
be expressed in a (d-l- l)-dimensional system by a Harris- 
type criterioni^, 

< 2/d . (1) 

Here iy± is the correlation length exponent in the spatial 
direction of the pure system . Indeed for directed perco- 
lation at any d < 4 dimensions the disorder is a relevant 
perturbation. 

The new, random fixed points, which control the criti- 
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cal behavior of absorbing state systems with quenched 
disorder have been studied in different papers and 
their properties are found in some respect to be 
unconventionali^ii^iiSiiLiSiiS. Early numerical studies by 
Noest^^ on random cellular automaton models in the di- 
rected percolation universality class in one and two di- 
mensions have shown considerable change of the criti- 
cal exponents in comparison with the pure system's val- 
ues. In more recent studiesi^ of the two-dimensional 
contact process with dilution, logarithmically slow dy- 
namical correlations were found. These could be re- 
lated to the results of a field-theoretical investigation 
by Jansseni^, according to which the renormalization 
group equations have only runaway solutions. Inter- 
estingly the exponents, associated with the logarithmic 
time-dependence of different dynamical quantities were 
found to be disorder dependent. The absorbing phase 
was found to have properties similar to that of a Grif- 
fiths phase and shows power law behavior with non- 
universal exponentsi^iiii^. The static critical behavior 
of the model has been explored less, but exponents in 
d — 2 have been determinedi^. 

In the present paper we revisit the problem of absorb- 
ing state phase transitions in the presence of quenched 
disorder. In this study we make use of the formal analogy 
between quantum systems and the " Hamiltonian" oper- 
ator formalism^SiSi of stochastic processes. Although the 
latter operators are generally non-hermitian, techniques 
developed in the study of quantum spin systems can often 
be successfully apphed^^-^^-^'^-^''. 

In the theory of random quantum spin systems, re- 
cently considerable progress has been made in under- 
standing their low-energy (or long time), long distance 
behavior. In particular, the use of a real space renor- 
malization group (RG) method, originally introduced by 
Ma, Dasgupta and Hu^^ has lead to many new, partially 
exact results and has given - at the same time - new 
physical insight into the problem^^. Subsequent analyt- 
ical and numerical work has provided further important 
results, in particular for one-dimensional system3SS*SL2£. 
One new concept which has emerged from these studies is 
the existence of strong disorder fixed points in which the 
disorder plays a completely dominant role. This property 
is manifested by the fact that during renormalization, the 
distribution of the random parameters (couplings, trans- 
verse fields, etc.) broaden without limits and therefore 
the RG treatment becomes asymptotically exact. 

Having the close similarity between quantum systems 
and the Hamiltonian description of stochastic processes, 
one might ask the question if the concept of strong dis- 
order fixed point applies for the latter in the presence 
of quenched disorder. For what is perhaps the simplest 
stochastic process, the random walk, the answer is posi- 
tive. Sinai diffusion^Si, which represents a particular type 
of random walk in a random environment, can be inter- 
preted as a realization of a strong disorder fixed point. 
Indeed RG methods have been successfully applied to 
explore new exact properties of this process^. 



While the Sinai diffusion obeys detailed balance, most 
of the stochastic particle systems do not. Therefore it 
is particularly interesting to know if quenched disorder 
could have a similar effect on the latter problems, too. 
In this paper we are going to study this issue in detail. 
In particular we investigate the applicability of the strong 
disorder RG scheme for models in different (non-random) 
universality classes. The RG predictions are then con- 
fronted with the results of extensive numerical calcula- 
tions, which are performed by density matrix renormal- 
ization (DMRG) and by Monte Carlo (MC) simulations. 
A short account of our results has been published in a 
Letter^. 

The structure of the paper is the following. Gen- 
eral notations about scaling theory at absorbing state 
phase transitions, both at conventional and strong dis- 
order fixed points, are given in section 2. The method 
of strong disorder renormalization group is explained in 
detail in section 3, where it is applied to the random con- 
tact process and to its generalization with two different 
absorbing states. Results on random directed percolation 
are obtained by mapping onto a random walk. Section 
4 contains a comparison of these analytical predictions 
with results of numerical calculations on the one- and 
two-dimensional contact process. Our conclusions are 
presented in the final section, and some technical details 
are given in the Appendices. 



II. SCALING AT ABSORBING STATE PHASE 
TRANSITIONS 

In the models that we consider here each site, i, of a 
d-dimensional lattice is either vacant (0) or occupied by 
at most one particle (A). The dynamics of the model is 
given by a continuous time Markov process and is there- 
fore defined in terms of transition rates. The reactions 
in the system are basically of two types: a) branching 
in which particles are created at empty sites (provided 
one of its neighbours is occupied) occurs with rate A^, 
and death of particles with a rate, /i^. The average num- 
ber of particles at site i, at a given time t is denoted 
by {ni){t). This system evolves to a stationary state in 
which averages become time independent. In this state, 
if the average value of the branching rates compared with 
the average value of the death rates is sufficiently large, 
a finite fraction of sites, p ^ ^ I L^^{ni) > 0, is occupied 
and the system is in the active phase. In the opposite 
situation, i.e. when the average value of the branching 
rates is relatively small, then p = and the system is in 
the inactive phase. The two phases described above are 
separated by a phase transition point, which is located 
at A = Ac, where A is a suitably defined control param- 
eter. In the homogeneous system with pi — p and A.; = A 
we have A = p/X. 
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A. Conventional scaling behavior 

In the vicinity of the nonequihbrium phase-transition 
point, where the reduced control parameter 6 = {A — 
Ac)/ Ac is small, generally anisotropic space-time scaling 
symmetry holds. The order-parameter, p{S, 1/t, as 
a function of a finite time-scale t, in a large finite system 
of size L, scales under changing lengths by a factor 6 > 1, 
V = L/h, 

p{5, lit, l/L) = b--p{6b'/'^- , 57t, b/L) . (2) 

Here the critical exponents, a:, and z have the usual 
definitions: the correlation length ^± in the spatial (the 
relaxation time tr or the correlation lenght ^|| in the tem- 
poral) direction diverges as ^ 1^1 {tr = ^|| 
|5|~'^ii) and the anisotropy or dynamical exponent is 
defined as z ^ v^^jv^. Taking the scaling parameter 
b = 6~'^-^ we have in the thermodynamic limit, l/L = 
and in the stationary state {1/t — 0), p{6) ^ 5^ with the 
order-parameter exponent (3 = xv^. For a surface site 
one defines the surface order-parameter, ps, the scaling 
behavior of which is the same as in equation Q but in- 
volves the surface scaling dimension, Xg, instead of the 
bulk exponent x so that Ps((5) ~ 5^= with (3^ = XsVat^- 

The dynamical behavior of stochastic systems is re- 
lated to the scaling properties of the characteristic time, 
tri which in the Hamiltonian formalism is given by the 
inverse of the smallest gap, e. In a conventional fixed 
point we have the relation: 

e{5, l/t, 1 /L) = b-^'eiSb^/"^ , b^t, b/L) . (3) 

thus in a finite system the appropriate scaling combina- 
tion is eL^. 

The autocorrelation function, G{5,l/t,l/L) = 
{ni{t)ni{0)) , which, at least for a homogeneous system, 
is independent of the position of a particle in the bulk 
has the scaling behavior at the transition point: 

G{6 = 0,1/t, l/L) = b-^^'GiS = 0, b'/t, b/L) . (4) 

Here, taking b — t^^^ we obtain in the thermodynamic 
limit G{t) ~ 

In a MC simulation one usually starts with one seed- 
particle in the origin of an empty lattice and measures 
the survival probability at the origin, Ps, the total num- 
ber of particles present in the system, N, and the mean- 
square distance of the particles from the origin, R^. For 
the contact process it can be shown^iS that in the long 
time limit, the survival probability equals the local order 
parameter. This 'duality' relation also holds in the disor- 
dered case. More generally, it is expected that if there is 
only one absorbing state, the two quantities are expected 
to obey the same scaling behaviour—. Thus from Eq.Q 
we then immediately get: 

Ps{S,l/t,l/L) = b--P,{6b'/''^,byt,b/L) , (5) 



For b = t^/^ we obtain at the critical point, Psif) ~ 
with 6 — x/z = P/{i/±z). Generally we have G{t) ~ 
Ps{t). The total number of particles is proportional to 
the integral of the density-density correlation function 
and obeys the scaling relation: 

N{5,l/t,l/L) ^ b'^-^^N{Sb^/''^,byt,b/L) , (6) 

It therefore behaves at the critical point as N{t) ~ t^, 
with Tj = [d— 2x)/ z. Finally, i?^ scales similarly as L^: 

R\S,l/t,l/L) = b^R\5b^/''^,byt,b/L) , (7) 

and at the critical point we have R^{t) ~ t^^^. We note 
that in the general situation, i.e when several absorbing 
states exist, the scaling relations in Eas. (|5l7(l involve a 
new exponent, x', instead of x. 

B. Strong disorder scaling 

Strong disorder fixed points were observed so far in 
random quantum systems as well as in Sinai diffusion. 
Here we suggest that the same type of fixed point can 
also exist in absorbing state phase transitions and we in- 
troduce the corresponding scaling theory, which is based 
on two ingredients^. First, at a strong disorder fixed 
point the dynamical behavior is ultraslow. The charac- 
teristic length-scale, ^, is related to the logarithm of the 
characteristic time-scale, t^, as: 

r^lnU . (8) 

thus the dynamical exponent, z, is formally infinity. Sec- 
ond, the average value of quantities related to the par- 
ticle occupation number are dominated by so called rare 
events (or rare regions of a large sample). In a rare event 
(ui) = 0(1), thus L independent, whereas in a typical re- 
alization (rii), is generally exponentially small in L. The 
fraction of rare events, /, is decreasing with the size, L. 
At the critical point we can write for it the scaling trans- 
formation: 

f{l/L) = b-^f{b/L) , (9) 

thus with b = L we have / ~ L^^. The average value 
of the particle density is indeed dominated by the rare 
events, thus p ^ f, and the order-parameter at the criti- 
cal point satisfies the scaling transformation: 

p{S,l/lnt,l/L) = b-'^piSb^/^^.b'^/lnt.b/L) . (10) 

One can see that the static exponents, (3 — xi^^ and 
f3s = x,sV}_ are given by the same expressions as in the 
case of ordinary scaling. (We shall explicitly show the 
construction of rare events for random directed percola- 
tion in section III.B2.) 

Due to the ultraslow dynamical behavior the scaling 
relation Q is modified into: 

lne(J,l/lni,l/i) = b'"^ \^e{ib^l''^ ,b^ /\^t,b/L) . 

(11) 
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thus in a finite system the appropriate scahng combina- 
tion is (Ine)L'^. 

To calculate the average autocorrelation function one 
should keep in mind that disorder in the time-direction 
is strictly correlated, thus in a rare event the autocor- 
relation function is of 0(1) and (almost) zero otherwise. 
Thus the average autocorrelation function is proportional 
to the fraction of rare events, G ~ /, and scales at the 
critical point: 

G{S = 0, 1/ \nt, l/L) = &-^G((5 = 0, b'*' / In t, h/L) . 

(12) 

Thus in the thermodynamic limit: Git) ^ (Ini) . 
We can also determine the scaling behavior of other dy- 
namical quantities such as Pg, N, and i?^, in the strong 
disorder fixed point in a similar way. It is therefore suf- 
ficient to replace in the conventional scaling relations in 
Eas (|5l7f) t by Int and z hy ip. As a consequence the 
time-dependence of the dynamical quantities will also be 
logarithmic at the critical point: 

Psit) ^ (Int)-' , 0^x/iP 
N{t) - (Int)^ , ri={d~2x)/^j 
R^(t) - (Int)^ , a^2/i) . (13) 

This type of logarithmic time dependence has been ob- 
served by Dickman and MoreirE^i^ by analyzing numerical 
data on the diluted 2d contact process, and were inter- 
preted as a "violation of scaling". The measured expo- 
nents 6, fj, and ct were found to be dilution, i.e. disorder, 
dependent. 

In the following, using a real space renormalization 
group method we i) shall give a natural explanation of the 
observed logarithmic time dependence and ii) shall calcu- 
late the critical exponents, which - in the one-dimensional 
case - are presumably exact. 

C. Scaling in the Griffiths-phase 

In a disordered non-equilibrium system, which is glob- 
ally in one stationary phase, say in the non-active phase 
with 5 > 0, there are specific local regions of size Zc, 
in which strong fluctuations of the local rates prefer 
the existence of the other phase, say the active phase. 
These rare regions, which are localized and have an ex- 
ponentionally small probability of occurance, p{lc) ~ 
exp(— a^c), contribute to an exponentially large relax- 
ation time^, tr ~ exp(tT/c). (In a d-dimensional sys- 
tem Ic should be replaced by l*^.) Then the distri- 
bution of large relaxation times has an algebraic tail: 

p{ir) ^ tr ^ , with 1/z' — 13/ a and the average auto- 
correlation function: 

G{t) - j dtrp{tr)exp{~t/tr) - t"^/''' (14) 

also decays algebraically. Consequently in this so called 
Griffiths phase^^ dynamical correlations are quasi-long- 
ranged, whereas spatial correlations are short ranged. 



The dynamical exponent is a continuously varying func- 
tion of the distance from the critical point, z' — z'{S). 

This result can be incorporated into a scaling theory 
as follows. Since a rare event, which brings the dominant 
contribution to the average autocorrelation function is 
localized, its probability of occurance is inversely propor- 
tional with the size of the system, L. Consequently the 
average autocorrelation function obeys the scaling law: 

G{5, l/t,l/L) = b-^G{S, b'^'/t, b/L) , (15) 

and with b — t^^^ we recover in the thermodynamic limit 
the relation in Ea. H14|l . Finally, scaling of the lowest gaps 
follow the rule: 

e{5, l/t,l/L) = 6-^'e(^, 6^' A, b/L) . (16) 

thus in a finite system the appropriate scaling combi- 
nation is eL^ . We note that in the Griffiths phase the 
power-law singularities are often supplemented by loga- 
rithmic correctiongSSi^, which are related to the fact that 
the size of the rare event grows logarithmically, Ic ~ InL, 
since p{lc) ~ l/L. 



III. THE RENORMALIZATION GROUP 
FRAMEWORK 

In this section we apply a real space renormalization 
group method for random stochastic particle systems, 
which is a variant of the Ma-Dasgupta-Hu method origi- 
nally developed to study random quantum spin chains^*. 
The essence of the method can be summarized in the 
following points. 

i) Start with the initial distribution of the random re- 
action rates, Pin(At) and i?in(A) and sort the rates in de- 
scending order. The largest rate, i.e. the fastest process, 
sets the energy scale, il, in the problem. 

ii) Integrate out the fastest local process, i.e. eliminate 
the rate, fl. This amounts to decimate out one site of the 
lattice or to replace a pair of sites with a new effective 
one. Renormalised couplings are then determined using a 
second order perturbation calculation. It is interesting to 
remark that for stochastic systems this step corresponds 
to a fast rate expansion as discussed in section 4.3 of 
reference 21. 

iii) Iterate the decimation process. This will result in 
a reduction of the energy-scale, fl, and a modification 
of the distribution of the (effective) rates: P(/i, 5^) and 
i?(A,rj). 

iv) At the fixed point of the transformation (which is 
a.t n = n* = 0) the distributions of the rates become sin- 
gular and from these singularities the value of the critical 
exponents are calculated. 

In the following we construct and solve the RG equa- 
tions explicitly for the random contact process. 
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A. The random contact process 

1. The Hamiltonian formalism 

In the contact process each site of the lattice can be 
either vacant (0) or occupied by at most one particle (A), 
and thus can be characterized by an Ising-spin variable, 
(7i = 1 for and (7i = —1 for A. The state of the system 
is then given by the vector P (cr, t) which gives the proba- 
bility that the system is in the state a = {. . . , a^, . . .} at 
time t. A particle can be created at an empty site i with 
a rate pXi/po, where p (po) is the number of occupied 
neighbors (the coordination number of the lattice) and 
at an occupied site the particle is annihilated with a rate 
fii. The time evolution is governed by a master equation, 
which can be written into the form: 



dP(^) 
dt 



= -HcpP{a) 



(17) 



Here the generator Hcp of the Markov process is given 
by: 



Hcp = ^ ^J.^Mi + ^ — {riiQj + QiUj) 



(18) 



in terms of the matrices: 



M 



-1 
1 



, n 




1 



,Q 



1 

-1 



and (ij) stands for nearest neighbors. It is well 
known^^ that the steady state probability distribution of 
a stochastic process coincides with the ground state of its 
generator (sometimes also called quantum Hamiltonian 
of the stochastic process) while relaxation properties can 
be determined from its low lying spectrum. 

For non-random couplings there is a non-equilibrium 
phase transition in the contact process which belongs to 
the universality class of directed percolationSi^. In one 
dimension it is at (/i/A)c = 0.3032 and the critical expo- 
nents are given by: /3 = 0.2765, Ps = 0.7337, v± = 1.097 
and z — 1.581. In two dimensions {ij./X)c = 0.6065 and 
[3 = 0.584, 13s = 1.03, va_ = 0.734 and z = 1.76. In the 
following we often use the variable A = A/po to charac- 
terize the creation rate. 



2. Decimation rules of the random process 

For the random contact process the transition rates 
and Ai are independent and identically distributed vari- 
ables and, as described in the previous section, the largest 
one Q. = max({Ai}, {/ii}) sets the energy scale in the sys- 
tem. Thus, the largest rate can be either one of the death 
rates /i^ or be one of the branching rates A^. For each a 
different way of decimation should be used. 

i) The largest term is a branching rate: — X2 



If the largest rate is a branching, say A2, which con- 
nects sites i = 2 and j = 3, the two-site cluster (2, 3) 
spends most of the time in the configurations AA or 00 
and can be rarely found in one of the other two configu- 
rations, A% and %A. Consequently for large times the two 
sites behave as a cluster with a moment of m = 2 and 
with an effective death rate, jl2, which can be calculated 
from the energy spectrum of the two-site Hamiltonian, 
H^p . As shown in the Appendix A the two lowest energy 
levels of the cluster are separated from the two others by 
a distance of A2, therefore in a good approximation just 
the two lowest levels can be retained. The effective death 
rate, jl, is given from the value of the lowest gap as: 



A* 



2/i2M3 

A2 



(19) 



where ^2 ^^nd /i3 are the original death rates at sites i — 2 
and J = 3, respectively. The renormalization equation 
in Ea. (|19|l should be extended by the renormalization of 
moments (i.e. the number of original sites in the cluster): 



m = m2 + , 



(20) 



where in the initial situation 7712 = ms — 1. 

The renormalized value of the death rate can also be 
obtained with the following reasoning. Let us start with 
the original representation, when the two-site cluster is 
in the occupied state, AA. In the effective decay pro- 
cess first the particle at site 2 should decay (with rate 
II2), which is then followed by the decay of the particle 
at 3. This second process has a very low probability of 
1^-3 /{X2 + Ma)- Since the same processes can also occur 
with the role of 2 and 3 interchanged, we find that for 
A2 ^ M2, the effective decay rate is given in Ea. (|19|l . 

ii) The largest term is a death rate: = ^2 
In this case the site (2) is almost always empty, 0, 
therefore it does not contribute to the fractal properties 
of the A cluster and can be decimated out. The effective 
branching rate, A between the remaining sites 1 and 3 can 
be obtained by studying the eigenvalues of the three-site 
Hamiltonian, H^p. As shown in the Appendix A out of 
the eight eigenvalues there are four of 0(/i2), which are 
discarded. The remaining lowest four levels are identified 
as the spectrum of a two-site cluster with an effective 
branching rate: 



A 



A2A 



2 A3 



M2 



(21) 



As in the previous case one can obtain the renormal- 
ized value of the branching rate with a simple reasoning, 
which works as follows. Let us have the configuration 
of the three-site cluster in the original representation as 
yl00. The effective branching rate between sites 1 and 3 
is generated by a virtual process, in which first a parti- 
cle is created at site 2 (rate A2), and then one at site 3 
(probability X3 / {X3 + 112)) ■ Hence, we get for very strong 
disorder the branching rate given in Ea. (|21|l . 
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The renormalization equations in Eas. H19|l and H21() 
can be transformed into a symmetric form in terms of 
the variable, J = A/k = A/ (po^) with k = \/2 as 



MM 



J = K— -. (22) 



We can see from Eg . (1221) that for weak disorder the gen- 
erated new rates can be occasionaUy larger then the deci- 
mated ones, thus in these steps the energy-scale does not 
lower. For strong enough disorder, however, these non- 
monotonic steps are expected to be so rare that they do 
not influence the behavior of the RG flow. With this as- 
sumption we analyze in the following the properties of 
the RG equations in one dimension. The results are then 
confronted with numerical calculations in section IV. 



3. Renormalization in one dimension 

In one dimension the topology of the lattice does not 
change under renormalization which makes it possible 

I 



to treat the problem analytically. First we note that 
after a repeated use of the transformations in Ea. l|22|) 
the generated branching (death) rates are in the form of 
a ratio of products of original branching (death) rates and 
original death (branching) rates. The control parameter, 
S, is defined as: 



[ln/x]av - [In J] av 
var[ln fi] + var[ln J] 



(23) 



and at the fixed point, 5 = 0, which follows from duality 
of the RG equations in Ea. l(^ (here [.]av denotes the 
average over disorder). 

At a given energy scale, we have the distribu- 
tion function of the death rates, P{n,ft), and that of 
the branching rates, R{J,fl). Changing the energy 
scale, ft n — dO, amounts to eliminate a fraction 
of drJ[P(n, n) + R{n, fl)] sites. The distribution of the 
branching rates changes as: 



R{j,n) + dnp{n,n) j ^Ji j dJ3R{Ji,n)R{J3M)^ 
{i-dn[p{n,n) + R{n,n)]y^ . 



S[J~^^]^S{J-Ji)-6{J-J3) 



(24) 



r 



Here in the r.h.s. the three delta functions represent the 
one generated new branching rate and the two decimated 
branching rates during one RG step and the second factor 
ensures normalization.'^^'.. A similar equation is obtained 
for the distribution of the death rates. From the duality 
of the RG equations, it follows that one should only make 
the interchange <-> J and P ^ R. 

Expanding R(J,n — dfl) one arrives to the integro- 
differential equation: 



^^R{j,n) [P{n,n)-R{n,n)] 

1^ ' 



p{n, n) I dj'R{j\ n)R [ -^n, n 



J/k 



J'k 



(25) 



and similarly one obtains for the distribution P(/i, fl): 
dP 



p{tj.,n) [R{n,n) - p{n,n)] 
-Rin, n) [ d^i'Pifi', n)p ( -^n, n] — . 



dfi 



(26) 



Solution of these equations can be obtained analytically 
at the fixed point, = 0, at 5 = 0, in which the dis- 
tributions i?( J, fi) and P{^^ f2) are asymptotically iden- 
tical. The calculations can be found in the Appendix B. 



According to these results the appropriate scaling vari- 
able is in logarithmic form rj = ~ (In — In J) / In il — 
— (Infi — ln/i)/lnil, and its distribution is given from 
Eqs.||BlJ and as 



p{r])dri = exp(— ?/)dr/ 



(27) 



The distribution in terms of the original variables, J (and 
/x) in Ea. HB2p is given by 



R{j,n)^R[^ 



1-Rn 



Rn = 



ln(17o/f))' 



(28) 



where fio is a reference energy scale, and the distribution 
becomes singular at the fixed point, as ^ 0. Due to 
this singularity the decimation transformation in Eq. H22() 
becomes exact at the fixed point. This can be shown by 
calculating the probability that one of the neighboring 
death rates, beside the largest branching rate with J = fl, 
has a value of /i > afi, with a < 1: 



P{a) ~ / P{J,a}dJ = Rn 
Jan 

« Rn\n{l/a) , 



-i+Rn 



dx 



(29) 
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which indeed goes to zero as the iteration proceeds, since 
— > 0. Consequently the RG transformation becomes 
asymptotically exact and the singularities, calculated by 
this method at the critical point, are also very probably 
exact. 

We start to determine the relation between the en- 
ergy scale, il, and the length scale, Lji, by studying the 
fraction of non-decimated sites, n^. When the energy 
scale is decreased by an amount of dfl a fraction of sites. 
driQ = nn[P(0) -I- R{^)], is decimated out, so that we 
obtain the differential equation: 



= nn[P{n) + R{n)] , 



which can be rewritten as 
din rifi 



dlnn 



n[p{n) + R{n)] = -2y{n) 



(30) 



(31) 



Using the solution to y(il) in Ea. ljB8|) one can integrate 
Eq.lPU with the result: 



nn = 



1 + 2/0 In^ 



-2 











.5 = 



(32) 



Thus from Eq. (|32|l we get for the typical distance between 
remaining spins, Lq, as: 



1 



In 



(5 = 



(33) 



This is equivalent with a logarithmic dynamical scaling 
as written in Eq.® with an exponent, i/' — 1/2. 

In order to calculate the singularity of other quantities, 
such as the correlation length and the order-parameter, 
at the fixed point we should study scaling of the lengths 
and the cluster moments. As we have shown in the Ap- 
pendix C, at the fixed point these calculations are equiv- 
alent to that for the random transverse Ising spin chain. 
Therefore here we quote only the results. For a detailed 
derivation we refer to the original literatur o^^i^^ . 

We remind first that a renormalized site is composed 
from parts of the original lattice and the renormalized 
length is given by the sum of the lengths in the original 
lattice. In the paramagnetic phase the average length of 
sites approach a finite value, during renormalization 
as f2 ^ 0. In the vicinity of the fixed point the RG- 
equations lead to a singularity: 
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(34) 



To study scaling of the order-parameter one should in- 
vestigate the average cluster moment, which at the crit- 
ical point behaves as: 



m 



mo 



In 



1 l + \/5 



(35) 



Then the order-parameter can be calculated as p = 
rn/Lfi, which behaves asymptotically as 



L- 



(36) 



and the scaling dimension, x, is given by Eqs.lO and 

as 



2-$ 



(37) 



Finally, scaling of the surface order-parameter is related 
to the average cluster moment of the surface site, mj. 
This is naturally smaller than for a bulk site, since the 
surface moment can accumulate from original sites only 
in one direction. According to an analysis of the RG 
results we have: 



1 I ^" 



thus 



1/2 



(38) 



(39) 



4. Renormalization in two dimensions 

In higher dimensions the topology of the lattice 
changes under renormalization: contacts and therefore 
reactions are generated between remote sites, too. How- 
ever, the renormalization does not introduce new types of 
reactions. Therefore the renormalization process, which 
is summarized in the decimation equations in Ea. l|22|) 
can be implemented numerically. As in one dimension 
for weak disorder the generated new couplings are fre- 
quently larger than the decimated ones, therefore the RG 
scheme does not work and scaling in the random system 
is most probably conventional, as described in section 
2A. For stronger disorder, however, the situation could 
change and strong disorder scaling could set in. If the 
critical behavior of the system is indeed attracted by a 
strong disorder fixed point then, as in one dimension, the 
value of a finite pre- factor k > in Eq. H22|) does not mat- 
ter. Numerical renormalization group calculations for the 
random transverse-field Ising model, for which k = 1, 
have shown the existence of a strong disorder fixed point 
in two-dimensions^Mi. The numerically observed crit- 
ical exponents are: x = 1.0, i'± = 1.07 and ip — .42 
Karevski et alm^ use a somewhat different numerical 
technique and find: x = 0.97, = 1.25 and V' = -5. 
In the light of the above arguments we expect also for 
the two-dimensional random contact process strong dis- 
order scaling with the above exponents in Ea. (|10|l if the 
strength of disorder is sufficiently large. To verify this 
scenario we shall reanalyze the numerical results of Mor- 
eira and Dickmani^ in section IV. B. 



B. Random directed percolation 

Directed percolation can be viewed as an anisotropic 
variant of percolation, in which the possible path of oc- 
cupied sites follows a given preferential direction. Equiv- 
alently, directed percolation can be interpreted as a 



8 



dynamical process, in which the spreading of a non- 
conserved agent is studied. In the random version of the 
problem the occupation probabilities are random vari- 
ables, which are strictly correlated along the preferential 
(time) direction, corresponding to random reaction rates 
in the dynamical interpretation. 

1. Random Reggeon field theory 

The field theory, which is expected to describe the 
critical behavior of directed percolation is the Reggeon 
field theory. In (1 -|- l)-dimension in the Hamiltonian 
limit the time-evolution of the process is governed by the 
Hamiltonian^ Hrf = J2i with: 

= -^.^f "I [(1 - 2a+)(l - 2<i) - af'Tf+i] (40) 

where cr^'^'^ are Pauli matrices at site, i, and — (erf -I- 
iaf)/2. The structure of Hup is similar to the Hamil- 
tonian of the random transverse-field Ising model in 
Ea. (|18|l : it consists of an interaction term (which is non- 
hermitian in Hjip) and a transverse field term. The order 
in the system is measured by the asymptotic limit of the 
autocorrelation function, G{t) = [(0|(Tf (t)crf (t)|0)]av: it 
is zero in the paramagnetic phase and finite in the or- 
dered phase. 

The renormalization of the random Reggeon-field the- 
ory can be made in a way that is completely similar 
to that for the random contact process (or the random 
transverse- field Ising model). A very strong coupling, 
f/2 = f^, will result in a two-site cluster in a renormalized 
field which is given by: 

(41) 

92 

On the other hand a site on which a very strong field, 
h2 = fl, acts is decimated out and a new coupling is 
generated between remaining sites, as: 



Comparing the decimation equations in Eqs.lQl} and H42() 
with those for the random contact process in Eas. (|19|l 
and H21(l we can see that they are equivalent. This is not 
surprising, since the two Hamilton operators in Eas. (|18|) 
and (|40|l are related through a unitary transformation'^, if 
yUi hi and Xi ^ gi- Consequently the singular behavior 
of the random contact process and the random Reggeon 
field theory are equivalent. 

Next, we consider the geometrical interpretation of di- 
rected percolation and study its critical behavior in the 
very strong disorder limit. 

2. Strong disorder: mapping to random walks 

Here we consider directed percolation on the square 
lattice with random occupation probabilities which are, 




q q q q q 



FIG. 1: Directed percolation on the square lattice with ran- 
dom bond occupation probabilities, which are perfectly cor- 
related along the vertical diagonal direction. Occupied bonds 
are drawn in bold. In the two columns on the left bonds are 
occupied with a large probability q, while in the other three 
columns this probability is small and equals q. The typical 
length-scale in the problem is I ~ 1/q. The direction of time, 
t, is also indicated which gives a reinterpretation of directed 
percolation in terms of a reaction-diffusion model. 



however, strictly correlated in the same layer, as shown 
in Fig. [T] 

For simplicity, we use a bimodal distribution, the oc- 
cupation probability in the i-th layer can be either pi = q 
with probability n or pi = 1 — q = q, with probability 
7f (we take q <q). For isotropic percolation the critical 
point is located at7r = 7f=l/2, which follows from self- 
duality^^. Since directed clusters generally contain less 
bonds than the isotropic ones, for directed percolation at 
the critical point tt < tt. In the strong disorder limit we 
take q —> in which case at the critical point {'k/tt)c — *■ 1 
as we shall show below. The effect of the same type of 
disorder for isotropic percolation has already been consid- 
ered by one of ua^. The reasonings in the two problems 
are very similar and they lead to identical critical proper- 
ties. Therefore, here we just briefly describe the method 
and emphasize the relations and differences between the 
directed and isotropic problems. For more details we re- 
fer to the original referencc^'^. 

The characteristic structure of the occupied bonds is 
very different for the two types of layers. In layers with a 
large probability, pi = q, almost all bonds are occupied, 
between two unoccupied bonds there is a characteristic 
distance, I ~ l/q. On the other hand in layers with a 
small probability, Pi = q, there are only very few occupied 
bonds, two occupied bonds are separated by the same 
characteristic distance, I ~ l/q. Note the duality: for 
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layers with pi = q the non-occupied (say white) bonds 
play the same role as the occupied (say black) bonds for 
layers with pi — q. 

Now let us consider a system consisting of r = 
1, 2, . . . , L layers with free boundary conditions, with the 
seed at r = and determine the surface order-parameter, 
Ps{L), which is given by the probability that the percolat- 
ing cluster extends up to the other surface at r = L. To 
calculate Ps{L) we consider strips of width k = 1,2, ... L 
and estimate, n{k), the typical number of sites in the 
k-the (i.e. surface) layer which are connected to the 
seed. It is evident from the duality properties that for 
n{k) ^ 1 the probability that the seed is connected to 
the k-th layer through white bonds is Pw(k) ~ l/n(k), 
since typically one out of n{k) sites has this property. 
We show by induction that n{k) is either zero or given 
by n{k) ^ q^^'-'^\ where X{k) > is the number of lay- 
ers with a probability q minus the number of layers with 
a probability q. First, our statement is trivially true for 
k — 1. Here n(l) ^ l/2g (and thuspiu(fc) ^ q), forpi = q 
and n{l) = 0, for pi ~ q. Evidently, if in a given sample 
n{k) — for some k < L, then n{k') = for any k' > k 
and the surface order-parameter is zero. To complete 
our proof, in the second step we show that, if X{k) > 1, 
then X{k + 1) = X{k) ±1, where the upper (lower) sign 
stands for a probability Pk+i = q [Pk+i = <?)■ The proof 
of this statement follows from the fact, that for Pk+i = q 
the number of black bonds in the cluster at layer k is 
reduced by a factor of ~ q, whereas for Pk+i = q the 
probability Pw{k) is reduced by a same factor ^ q. Here, 
when X{k) — 1 and Pk+i — q, thus X{k + 1) = the 
cluster is considered to be terminated at this point, thus 
Ps{L) = 0. 

To calculate the average value of the surface order- 
parameter we use a random walk picture, which has al- 
ready been applied in the isotropic problem^. To a 
given sample with a given probability distribution we 
assign a random walk which starts at position Xq — 
and makes its i-th step upwards (downwards) for pi — q 
(pi — q). The position of the walker in the k-th step is 
just Xk. The existence of finite surface order-parameter 
in the given sample is then formulated by the condi- 
tion: n{k) > for /c = 1,2, . . . ,L, thus Xk > for 
k = 1,2, . . . , L so that the random walk has a surviv- 
ing character. The average value of the surface order- 
parameter is given by the fraction of samples with finite 
surface order, which is just the survival probability of 
the random walk: [ps]av{L) ~ Ps{L) ^ L~^/^ . From this 
relation the value of the surface order-parameter scaling 
dimension, Xs = 1/2 follows, which is the same as for 
the random contact process in the strong disorder fixed 
point in Eq. (|^ . Note that the rare events introduced in 
section H.B in this case are the samples with a surviving 
characteristics in the probability distribution. Further- 
more in the critical situation, when the average surface 
order-parameter decays as a power with L the probabil- 
ities, TT = tF, as announced before. 

Other exponents can be deduced in analogous way as 



for isotropic percolation. The typical temporal extent 
of the percolating cluster is given by the typical num- 
ber of connected sites in a layer, ~ ntyp{L) ~ g"^'"", 
where the typical excursion of a (surviving) random walk 
in L steps is Xtyp ^ L^^^. Thus we have the logarith- 
mic scaling relation in Eq.JHJ with il) = 1/2, as for the 
random contact process. To calculate the bulk order- 
parameter the seed is put in the middle of the lattice and 
in a given sample p{L) — 0{1) if the percolating cluster 
has an extent of L. In the language of random walks this 
property is related to the so called average persistence^^. 
From this, scaling of the average (bulk) order-parameter 
IS given by [p]av(i) ~ L-"", with x = {3 - V5)/4, as for 
the random contact process in Eg. 13711 . Finally, to de- 
termine the scaling exponent of the order-parameter in 
the vicinity of the critical point one should consider ran- 
dom walks with finite bias, (5^ 7^ 0, towards the absorbing 
walk and S S^. From the scaling of the surviving prob- 
ability of biased walks'^'' one obtains for the correlation 
length, ^ I'^l^^- Thus we recover the same exponent, 
i'± = 2, as for the random contact process in Ea. (|34|) . 

Hence we can conclude that the critical properties of 
random directed percolation can be deduced from a ran- 
dom walk mapping in the strong disorder limit. The 
critical behavior is the same as for random isotropic per- 
colation and corresponds to the RG results obtained for 
the random contact process in the strong disorder fixed 
point. 



C. The generalized contact process with disorder 

So far we have considered different variants of absorb- 
ing state phase transitions which, in the absence of dis- 
order, all belong to the directed percolation universality 
class. We observed, that in the presence of strong enough 
quenched disorder all these processes show strong disor- 
der scaling behavior with identical critical exponents. In 
this section we consider other processes, which are not in 
the directed percolation universality class. 

For a renormalization group treatment, one of the 
most convenient models is the generalized contact pro- 
cess with several different absorbing states introduced 
by Hinrichsen^. In this model a site can be occupied 
by a particle. A, or can be in one of n empty states 
0i,02,...0n. Furthermore, besides the rules known for 
the ordinary contact process, there is a competition be- 
tween the different type of empty states. At the border 
of clusters with different type of empty states particles 
can be created. As a consequence, in these models, with 
increasing n there is a preference for the active phase 
and the phase transitions is found to be in universality 
classes different from the directed percolation one. For 
n — 2, the model was shown to be in the parity con- 
serving class^, whereas for n > 3 the model is always 
active^. 

Here we consider in particular the effect of disorder 
on the model with n = 2. For this case, the following 
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processes are allowed: 

AA^ A%i,A%2,^iA,%2A rate ^l,/2 

^0i,0iA-> 0i0i; A02,02^^ 0202 rate ^ii 

^0i,A02,0iA02^-> rate A, 

0102-^01^,^02; 0201 ^02^,^01 rate A, (43) 

For non-random couplings the phase-transition in one di- 
mension is located at /i/A = 1.592 and the critical expo- 
nents are = 1.82, z = 1.75, /3 = 0.91 consistent with 
those of the parity conserving universality class. Note 
that the Harris-type criterion in Eq.l^ with the above 
i^_L predicts a relevant perturbation for quenched disor- 
der. 

Next, for the random system we try to apply the renor- 
malization group method along the lines used for the ran- 
dom contact process in section III. A. We start with the 
decimation scheme and consider the situation when the 
largest rate is a branching rate, say A2 = As for the 
random contact process we take the two-site Hamilto- 
nian, which contains now 9 states, since each site could 
be in three (one active and two different inactive) states, 
and calculate the lowest three eigenstates along the lines 
presented in the Appendix A for the random contact pro- 
cess. Among the three lowest states, which are kept af- 
ter decimation, two have eigenvalues, 0, and the third 
has an energy ef' ~ 2/^2 /is /A2. Thus after decimation 
we have an effective two-site cluster in the presence of a 
renormalized death rate, which is given just in the same 
form as in the random contact process, see Ea. H19l) . The 
value of the effective death rate can be also obtained by 
a similar argument as for the random contact process. 
In order to calculate the decay from AA — > 0i0i,0202 
one should consider two second-order virtual processes: 
AA A%i -> 0101 and AA A%2 0202, which leads 
to Ea. (|19(l by taking into account the definition of the 
rates in Eq. (031 • 

For the case of a strong elimination rate, ^2 — the 
three site Hamiltonian contains 27 eigenstates, which are 
divided into 9 orthogonal sectors, each of which have 3 
states. The highest levels of each sector have an energy of 
0(^2) and thus can be discarded during decimation. The 
remaining 18 states, however, are all in the same order of 
magnitude. Three sectors have ground state energy zero 
and first excitation energy, A2 -|- A3. Two other sectors 
have the lowest energies: 



,123 



[7A2 + 4A3) ± \/(7A2 + 4A3)2 - I6A2A3 



(44) 

and in another two we should exchange in Eg. H44|l A2 and 
A3. Finally the last two sectors are also degenerate with 
the lowest eigenvalues: 



,123 



3(A2 + A3) ± \/9(A2 + A3)2 - 32A2A3I /4 . (45) 



assigned to one renormalized site and ii) the remaining 
energy levels are of the same order as the original rates, 
thus the energy scale is not lowered during these steps. 
These features of the decimation can be seen by the fol- 
lowing argument, too. With a large death rate, ij,2 = O, 
the site i = 2 is almost always inactive, so that it is ei- 
ther in 01 or 02 most of the time. Suppose we decimate 
this site and calculate the effective rate for a process in 
which in the original lattice at sites 1 and 3 there is A 
and 01, respectively, which will change to A and A. In 
the original bases this process can most easily be realized 
through A|02|0i — > A|02|A, which i) has a rate A2 and ii) 
the effective rate does not depend on the fact that site 
1 is occupied. Consequently during renormalization new 
degrees of freedom will appear and there is no systematic 
decrease of the energy-scale. 

Thus we can conclude that the renormalization group 
scheme does not work for the generalized contact pro- 
cess. Therefore, it is improbable that a strong disorder 
fixed point with the scaling properties described in sec- 
tion II. B could be present in this model for some finite 
value of the disorder. The qualitative difference between 
the random contact process and the generalized random 
contact process is due to the competition between the dif- 
ferent absorbing states in the latter model. For this case 
a large death rate indirectly promotes particle branching, 
since the active phase intrudes between the different ab- 
sorbing states. This is the reason why the pure model is 
always in the active phase for n > 3 — . This could also 
be true for the random model. 



IV. NUMERICAL INVESTIGATIONS 

Here, we discuss the results of two numerical ap- 
proaches which allow us to investigate the disordered 
contact process as a function of disorder strength. In 
this way, we will be able to show that the predictions 
made in the previous section are consistent with the nu- 
merical results provided the disorder is strong enough. 
However, we also find that there is a regime at small to 
intermediate disorder in which the critical exponents de- 
viate from the strong disorder ones and indeed seem to 
vary continuously. This is evidence for a line of fixed 
points. However, the numerical results available do not 
allow us to decide whether these are ordinary disorder 
fixed points (finite z) or strong disorder ones (infinite z). 



A. One-dimensional random contact process 

In our numerical work, we investigated the particular 
case for which the rate /i^ = 1 and the branching rate A^ 
is distributed according to 



Consequently the decimation does not work out here for 
a large death rate, since i) one can only discard 9 out of 
the 27 cell states, so that the remaining states can not be 



i?(A) ^ [^(A-A+) + 5(A-A_)]/2 



(46) 



with A± = exp (A±VDy The main advantage of this 
distribution is that it allows an exact average over dis- 
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order to be made, at least on lattices that are not too 
large. For this choice of R{X) the average value of In A 
equals A whereas the variance is D. These can there- 
fore be considered as suitable parameters to measure the 
activity and the disorder respectively. 



1. DMRG studies 

The density matrix renormalisation group (DMRG) 
method is a numerical technique that was originally in- 
troduced to investigate the properties of quantum spin 
or fermion systems. The method allows a precise deter- 
mination of the properties (energy, magnetisation pro- 
files, correlations, ...) of the ground state and low lying 
excitations of such models. The method is most suc- 
cesfull in one dimension. Given the formal similarity 
between quantum systems and stochastic ones, several 
groups started to apply this technique to interacting par- 
ticle systems in recent ye&rs^^^. The method is now 
known to work well also in these cases though it cannot 
give as accurate results as for the spin chains, mainly 
because at this moment algorithms to diagonalise non- 
hcrmitean matrices are not as well developped as those 
for the hermitean case. 

In our DMRG work we made calculations for systems 
with up to L = 24 sites. For L < 14 we were able to 
perform an exact average over all possible realisations of 
the disorder. For the larger systems sizes we considered 
typically around 10* disorder realisations. These calcu- 
lations were possible for D < 2. For larger D-values, 
we encountered numerical difficulties in the DMRG algo- 
rithm. 

In any finite system, the stationary state of the contact 
process is the absorbing state, i.e. the lattice without any 
particles. In order to study the absorbing state phase 
transition within finite systems, we therefore worked with 
open boundaries and put the death rate at the most right 
site of the lattice, /i^, equal to zero. This ensures the 
presence of particles in the stationary state. It can be 
expected that for large enough systems and for sites that 
are deep in the bulk the effect of this boundary will be 
negfigible. 

Using the DMRG we then firstly calculated the ground 
state of the generator (|18|) with open boundaries. Within 
this ground state we calculated the density profile, i.e. 

Pi = [(ni)]av 

In order to determine the location of the critical point 
and the critical exponents x, Xg and i^j^ we investigated 
in detail the behaviour of ps = pi (surface density) and 
Pl /2 (which we took as our estimate for the bulk density) 
as a function of the parameters A en D. 

A major problem of the contact process in comparison 
with several well studied disordered quantum spin chains 
is that in the present case the location of the critical point 
is not known exactly. The numerical inaccuracy in the 
location of the critical point will in its turn inffuence the 
accuracy by which critical exponents can be determined. 
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FIG. 2: Scaling plot of the particle density at D = .25 assum- 
ing Ac = 1.19, X — .24 and u± — 1.25. The different symbols 
indicate L = 12(D), L = 14(0),^- = 16(A), L = 18{s^),L = 
20(0), L = 22(0) and L = 24([>). 



To locate the critical point we investigated the quantity 



Y = dL 

^ d\nL 



(47) 



For a homogeneous system, p reaches its large L-value 
exponentially fast away from the critical point, and as 
a power law at the critical point (see |(5J) and (jSjl which 
are the same in the stationary state). As a consequence 
Yl goes to —00 for a non-critical system and to —(1 -I- a;) 
at criticality. We can therefore expect that in a finite 
system, the quantity Yl goes through a maximum as a 
function of fi/X from which one can obtain a finite size 
estimate of the location of the critical point and of a;. In 
this way we have obtained the critical value of A, denoted 
as Ac for different values of the disorder strength D and 
for different L-values. An extrapolation for L —^ 00 then 
gives our final estimates for Ac{D). The data for the 
surface density ps can be analysed in a completely similar 
way and give an independent estimate for the location 
of the critical point. From the analysis in section II, we 
expect that for _D — > 00, the critical point obeys [In /i]av — 
[In J]av which for the present case leads to the prediction 
of the exact location of the critical point, at Ac = In a/8 = 
1.0397. Unfortunately, we are not able to investigate very 
large Z?- values and hence could not verify this prediction. 

To determine the stationary state critical exponents we 
next made scaling plots for the density and the surface 
density assuming ((SJ (taking t 00). In Fig. |2]we show 
such a scaling plot for the density and for D = .25 and 
Ac = 1.19. From this we obtain x = .24 and = 1.25. 
Data for the bulk and surface density for other values 
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of D can be analysed in a completely similar way and 
they give rise to the DMRG estimates for the exponents 
Xs and X as shown in Fig. 3 of reference 31 and Fig. El 
respectively. From these we see that Xs assumes the value 
predicted for the strong disorder fixed point at I? ~ 1.5. 
For smaller _D-values, Xs decreases continuously starting 
from its value for a homogeneous contact process {x^ ~ 
.669 iS). 

A rather similar behaviour is found for the bulk expo- 
nent X. For increasing D, its value decreases from that 
for the homogenous contact process, x ~ .252, to a value 
of X .2 in the region 1.5 < Z? < 2. Here it has to be 
remarked that since the DMRG can only be performed 
for systems with L up to 24, it may very well be that our 
estimate of the bulk density p is still influenced by the 
boundary conditions which we had to choose, and this 
may very well be the reason why the exponent x reaches 
its value at the strong disorder fixed point more slowly. 

The exponent is even more difficult to estimate 
since its value is very sensitive to the precise location of 
the critical point. The values that we have determined 
are only rough estimates. We find that the value of i^j. 
increases from the pure system value with increasing dis- 
order, and equals 1.67 ± .08 at D = 1.5. 

The exponent z, or in case of a logaritmic scaling ip, 
can in principle be determined from the distribution of 
the gap in the spectrum of the generator H18|l. Using the 
DMRG, we therefore also calculated the distribution of 
gap sizes for different system sizes. In Fig. O we show 
the results of such an analysis at the critical point and 
for _D = .5. In the upper graph we show a plot assuming 
ordinary scaling and the value z = 2., while in the lower 
part we have assumed logarithmic scaling and ip = .35. 
As can be seen from this figure, the DMRG-results do 
not allow us to discriminate between the two types of 
dynamical scaling. A similar conclusion holds for other 
-D- values. 



2. Monte Carlo simulations 

In an attempt to reach bigger system sizes and to get 
independent estimates for the exponents we also per- 
formed extensive numerical simulations taking a single 
seed particle as initial condition. In order to have a good 
comparison with DMRG-data, we again put all death 
rates equal to one, while the branching rates are dis- 
tributed according to • For various values of D and 
A we calculated the survival probability Ps{t), the total 
number of particles in the system N{t) as well as their av- 
erage spread R'^{t). We typically simulated up to t = 10^ 
and took an average over 2.10^ disorder realisations. For 
large -D-values, it is difficult to obtain reliable simulation 
data since the dynamics becomes extremely slow. For 
this reason we are not able to explore values of D that 
are much larger then 1.5. 

To analyse our data we first have to determine the loca- 
tion of the critical point Ac{D). While in homogeneous 
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FIG. 3: Test of dynamical scaling form for D = .5, Ac = 1.19. 
In the upper graph, we have assumed ordinary dynamical scal- 
ing with z = 2.0, while in the lower graph we used logarithmic 
scaling with ip = .35. The different symbols correspond with 
L = 18(v), L = 20(A), L = 22(0) and L = 24(n). 



systems the critical point is the only one characterised 
by a power law decay of Psit), the same is not true in a 
disordered system as first pointed out by Bramson, Dur- 
rett and Schonmannpi. Power law behaviour can indeed 
be found in the whole subcritical regime and is a man- 
ifestation of the presence of a Griffiths phase. Still, we 
expect N{t) to increase above Ac and to decrease below 
criticality. In this way we obtain a first estimate of the 
location of the critical point. A second criterion which 
we used is that both for convential and for strong disor- 
der scaling, at criticality, InP/lnN becomes a constant 
asymptotically. In Fig. 0] we show some typical set of 
data for this quantity taken at D = 0.5. From these we 
determine Ac{0.5) ~ 1.177. The values for the critical 
activities which we find in this way are close (approxi- 
mately within 1%) to these found from the DMRG. 

In Fig. |Slwe present log- log plots for the three quan- 
tities of interest at the critical point with D = 1. These 
results are typical also for the other D-values investi- 
gated. 

As can be seen there is still some curvature visible in 
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FIG. 4: Plots of In P{t) / In N{t) as a function of In t at D = .5 
and for A = 1.179, 1.177 and 1.175 (top to bottom). From 
this,we estimate Ac — 1.177. 
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FIG. 5: Plot of the survival probability Ps{t), the number of 
particles N{t) and their average spread R^{t) (bottom to top) 
as a function of t for D = 1, Ac = 1.171. 



FIG. 6: Plot of the exponent x versus disorder strength 
D. The different symbols indicate the results obtained from 
the DMRG (squares) , and the simulations assuming ordinary 
scaling (circles) and logarithmic scaling (triangles). We also 
indicate the values for strong disorder (dotted line) and for 
the homogeneous system (dashed line). 



these figures, which could e.g. arise from the logarithmic 
corrections which are ubiquitously present in these kind 
of random 'quantum' systems. Yet, the late time part 
can be fitted quite well to a power law from which we 
can determine estimates of the exponents x and z as a 
function of D. Our results for these exponents are pre- 
sented in Fig. Eland Fig. [71 respectively. 

We notice two interesting trends. Firstly, we observe a 
decrease of x from its homogenous system value to a value 
that is consistent with that at the infinite disorder fixed 
point. The numerical values are moreover consistent with 
those found from the DMRG. Secondly, we observe that 
the dynamical exponent z seems to make a jump as soon 
as any disorder is present after which it increases with 
D. These results are then consistent with the idea of a 
line of ordinary (disorder) critical points which ends at 
some critical value of the disorder above which exponents 
assume precisely the values of the strong disorder fixed 
point. 

In a regime governed by strong disorder fixed points 
we expect that the dynamical quantities do not follow a 
power law but scale logarithmically, i.e. as given in H13|l . 
Indeed such a behaviour was found in the simulations in 
d = 2 — . We found that this kind of scaling can also de- 
scribe our results for all Z?- values investigated. Assuming 
strong disorder scaling, we can determine estimates for 
the exponents x and ip. These values are also shown in 
Fig. ini and Fig. [3 We notice that the values for the 
exponent x are not very sensitive to the type of scaling 
that we assumed. We also observe that the exponent ip 
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FIG. 7: Plot of the exponents l/2:(squares) and ?/>(circles) as 
a function of disorder strength D. 



seems to increase towards its value expected at the strong 
disorder fixed point, i.e. ^ — 1/2. 

We can thus see that, as was for the case for the 
DMRG, it is not possible to rule out from the numerics 
that there is a line of strong disorder fixed points in the 
model. When we compare the estimates for the dynami- 
cal exponents coming from the two numerical approaches, 
we find that those for the logarithmic corrections are 
more selfconsistent. There is a rather large discrepancy 
between the z-values as found from the DMRG and the 
simulations. This could be due to a lack of asymptotic- 
ity in time for the simulations and in size for the DMRG. 
Assuming logarithmic scaling, the values which we find 
from the two approaches are almost the same. The ex- 
istence of a line of strong disorder fixed points was not 
found so far in real quantum systems and could be an es- 
sential new feature of stochastic models. Yet, at present 
we have no theoretical underpinning for the existence of 
such a line of disorder fixed points. Further numerical 
studies are needed to obtain more insight on this point. 



FIG. 8: Values for the exponents a;(squares) and ^(circles) 
as a function of dilution p assuming logarithmic scaling. The 
values in this table are for d = 2 and are determined from the 
numbers given in table I of the first paper in reference 15. 



dataset, it is not possible to investigate whether the data 
of these authors are also consistent with ordinary scaling 
in the small p-regime. From the numerically determined 
values of 9, rj and ct it is however possible to determine 
values for the exponents x and in the two-dimensional 
case. The results are shown in Fig. |S1 

These numbers can now be compared with those ex- 
pected to hold at the strong disorder fixed point of the 
RTIM in two dimensions, which are x = l.Q and %p = .42. 
Note that for the largest disorder value, the exponent ip 
is very close to this value. The exponent x seems larger 
then expected. However, in this respect it is interesting 
to remark that the authors of ref. 15 notice that the 
data on N{t) are consistent with ordinary scaling but 
with = or with logarithmic scaling and a very small 
value of rj. This is consistent with our ideas which pre- 
dict that at the strong disorder fixed point in d = 2 and 
for X = 1, rj equals zero. Hence, we believe that also in 
two dimensions the strongly disordered contact process 
is in the same universality class as the RTIM. 



B. Two-dimensional random contact process 

As already remarked above, in the second paper of 
reference 15, the logarithmic scaling (|13|l was first ob- 
served in the two-dimensional contact process with dilu- 
tion where it was dubbed 'violation of scaling'. In the 
present work we have shown that this kind of scaling 
is the natural one to be expected at a strong disorder 
fixed point. In reference 15 it was assumed that the log- 
arithmic scaling holds for all values of the dilution, which 
here we will denote by p. Since we don't have the original 



V. CONCLUSIONS 

In this paper, we have investigated the effect of 
quenched disorder in the transition rates on the criti- 
cal behaviour of some models with absorbing state phase 
transitions. For models in the directed percolation uni- 
versality class, we have given convincing evidence that, 
if the disorder is large enough, the universal behaviour is 
that known from the RTIM and some other disordered 
quantum spin chains. This result follows from calcula- 
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tions using a strong disorder renormalisation approach 
and is consistent with numerical results in one dimension 
obtained with the DMRG and simulations. A reanalysis 
of data for a diluted contacted process in two dimensions 
is also consistent with this conclusion. It therefore seems 
that the effects of strong disorder lead to a new kind of 
universality between hermitean and non-hermitean quan- 
tum models. 

For weak and intermediate disorder our numerical re- 
sults indicate the existence of a line of fixed points with 
continuously varying exponents. Such behaviour has 
been found previously in some quantum spin chains'*^. 
Most interestingly, our data could be consistent with a 
scenario in which this is a line of strong disorder fixed 
points. Such a behaviour is also consistent with data 
obtained by Moreira and Dickman on a two-dimensional 
diluted contact process. Further simulations are needed 
to confirm this picture. It would be especially interesting 
to develop techniques that can efficiently simulate these 
kind of systems for sufficiently large disorder and long 
enough times. Moreover, it would be very interesting to 
investigate whether the behaviour that we find can also 
be seen for distributions of the transition rates different 
from the ones that we used in our simulations. If the log- 
arithmic dynamical behaviour can be demonstrated, it 
becomes a challenge to understand such behaviour from 
analytical or RG arguments. 

For absorbing state phase transitions in other univer- 
sality classes we have fewer results. Our RG calculations 
indicate that for the generalised contact process with dis- 
order, it may be possible that the model is always ac- 
tive, also for n — 2. This prediction should be verifiable 
with simulations. Whether this result also holds for other 
models in the same universality class is not clear yet. We 
believe that the effect of quenched disorder on absorbing 



state phase transitions in particular, and on stochastic 
many-particle systems in general, provides an interesting 
field of future research. 



APPENDIX A: DECIMATION FOR THE 
RANDOM CONTACT PROCESS 

When the branching rate between two sites (say 2 and 
3) is the largest rate in the system, we consider the Hamil- 
tonian H18|l restricted to these two sites and with free 
boundary conditions. This operator can be represented 
by the matrix: 



-^3 


-/i2 





^3 + A2 





-M2 





("2 + A2 




-A2 


-A2 


/i2 + M3 



(Al) 



In the limit ^2/^2 ^ and IJ-3/X2 the eigenval- 
ues are 0, 0, A2, A2. The degeneracies disappear for finite 
/i2/A and ^i^/X. Discarding the two highest energy levels, 
the gap between the remaining two lowest states can be 
calculated perturbatively leading to: 

E'^P = (A2) 
A2 

Now keeping in mind that a one-site cluster with a death 
rate // has a gap /i we obtain for the renormalized death 
rate, fl2 = E^p, as announced in H19|l . 

When the largest rate is a death rate, say fi2, one con- 
siders a three site cluster whose Hamiltonian (free bound- 
ary conditions) is represented by the matrix: 





-IJ'2 



















/i2- 


f A2 + A3 


























A2 



















-A2 


-A2 


M2 + A3 


























A3 


-^J.2 











-A3 








-A3 


fJ-2 + A2 


























A2 + A3 


-Ai2 










-A3 





-A2 


— A2 — A3 


fJ'2 / 



(A3) 



The corresponding eigenvalue problem is reduced to the 
diagonalization of four 2x2 matrices. In each subspace 
there is an eigenvalue of 0(/i2), which is discarded. Of 
the remaining four lowest eigenvalues there are two zero 
and two with the value 

E'S = ^ ■ (A4) 



five branching rate A the spectrum is given by 0, 0, A, A 
we obtain for the renormalized rate A = E^p, the result 
stated in (|^ . 



If we keep in mind that in a two-site cluster with an effec- 
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APPENDIX B: SOLUTION OF THE RG 
EQUATIONS AT THE CRITICAL POINT 

We look for a special solution of (|25|l and (|26|l (in the 
fixed point 17 ^ 0) of the form: 



1-P(0,0)0 



Substituting (jBl|) and (|B2|) into we get as: 



(Bl) 
(B2) 



n — \ = P-R-QRP n k""^ , 

an dn J n j ' 

(B3) 

where we used the notations, = P{^, ^) and 

i?(ri) = i?(r2,r2). For a moment we set k = 1. After 
a trivial rearrangement, we obtain the relation: 



rtRP- 



dflR 

dn 



^ J nil 



which leads to the ordinary differential equations: 



dR 

dn 

dP _ 
dn 



R 

n 



PR 



p 



(B4) 

(B5) 
(B6) 



These equations, which hold just for k = 1 can be 
solved for general (non-symmetric) distributions^. The 
solution gives information about the singular behavior of 
the system outside the critical point, i.e. in the Griffiths 
phase. This general solution, however, does not apply for 
K 7^ 1. In contrary, the solution at the fixed point, when 
R{J,il) and P(/^,J7) are identical, thus R = P, holds 
for any finite k > 0. Indeed, in terms of the variable, 
y = Rn = Pil, and the log-energy scale, F = — In we 
obtain the differential equation: 



dy_ 
dF 



with solution: 

y = Rn = Pn 



F - Fo ln(r2o/f^) 



(B7) 



, ^ = , (B8) 



Here Fo = — Infio is a reference (log)energy scale. Now, 
going back to (|B3p we can see, that at the fixed point the 

actual value of k does not matter. The term k~^^ 1 
and \nK/lii{n/J) ~ kM 0. 



APPENDIX C: THE RANDOM 
TRANSVERSE-FIELD ISING MODEL AND ITS 
DECIMATION RULES 

The prototype of random quantum systems is the ran- 
dom transverse-field Ising model which is defined by the 
Hamiltonian: 



H 



(CI) 



Here the sum runs over nearest neighbors and erf, af 
are Pauli matrices at site i. The exchange couplings 
Jij and the transverse-fields hi are independently dis- 
tributed random variables with distributions 7r(J) and 
p(h), respectively. The Hamiltonian in Ea. ljCl|l in one 
dimension is closely related to the transfer matrix of a 
classical two-dimensional layered Ising model, which was 
first introduced and studied by McCoy and Wu^. 

In the renormalization scheme we have the following 
decimation relations. If the largest term in the Hamil- 
tonian is a coupling, say J — Q, then the two sites con- 
nected by J and having transverse fields, h and h' and 
moments m and m', behave as an effective composite 
cluster with moment rn in a renormalized transverse field, 
h, which are given by: 



h = 



hh!_ 



rh = m + m 



(C2) 



On the other hand if the largest term is a transverse 
field, say /i = fi, then the site with this transverse field 
gives negligible contribution to the (longitudinal) mag- 
netic susceptibility, thus can be decimated out. This 
leads to a renormalised coupling between the nearest 
neighbours of the decimated site that is given by: 



J 



JJ' 



rh = m + m 



(C3) 



which is related to 1C2|I through duality. 
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